This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
library(spdep)
## Loading required package: sp
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.5.1, GDAL 2.2.2, PROJ 4.9.2
library(maptools)
## Checking rgeos availability: TRUE
library(leaflet)
chi.poly <- readShapePoly('./shapefiles/PovertyData.shp')
## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read
class(chi.poly)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
str(slot(chi.poly,"data"))
## 'data.frame': 88 obs. of 27 variables:
## $ OBJECTID_1: int 1 2 3 4 5 6 7 8 9 10 ...
## $ COUNTY_CD : Factor w/ 88 levels "ADA","ALL","ASD",..: 36 37 38 39 25 40 41 42 61 63 ...
## $ COUNTY_SEA: Factor w/ 88 levels "AKRON","ASHLAND",..: 33 42 51 58 20 35 70 53 10 61 ...
## $ ODOT_DISTR: int 9 10 11 3 6 9 11 5 10 1 ...
## $ FIPS_COUNT: Factor w/ 88 levels "39001","39003",..: 36 37 38 39 25 40 41 42 61 63 ...
## $ POP_2010 : int 43589 29380 42366 59626 1163414 33225 69709 60921 14645 19614 ...
## $ POP_2000 : int 40875 28241 38943 59487 1068978 32641 73894 54500 14058 20293 ...
## $ POP_1990 : int 35728 25533 32849 56240 961437 30230 80298 47473 11336 20488 ...
## $ STATE_PLAN: Factor w/ 2 levels "N","S": 2 2 1 1 2 2 1 1 2 1 ...
## $ ELEVATION_: int 1340 1220 1380 1200 1130 1040 1420 1420 1340 780 ...
## $ ELEVATION1: int 700 640 800 600 670 600 640 840 660 650 ...
## $ LAT_NORTH_: num 39.4 39.7 40.7 41.3 40.1 ...
## $ LAT_SOUTH_: num 39 39.4 40.4 41 39.8 ...
## $ LONG_EAST_: num -83.3 -82.2 -81.6 -82.3 -82.8 ...
## $ LONG_WEST_: num -83.9 -82.7 -82.2 -82.8 -83.3 ...
## $ AREA_SQMI : num 558 424 424 496 544 ...
## $ AREA_ID : int 53 54 55 56 143 57 58 59 75 76 ...
## $ SHAPE_STAr: num 2.41e+09 1.84e+09 1.91e+09 2.27e+09 2.40e+09 ...
## $ SHAPE_STLe: num 213158 215452 193045 197349 211754 ...
## $ NoOfUni : int 0 0 0 0 13 0 1 2 0 0 ...
## $ Shape_Leng: num 1.72 1.75 1.59 1.58 1.7 ...
## $ Shape_Area: num 0.151 0.115 0.117 0.138 0.149 ...
## $ FIPS : int 39071 39073 39075 39077 39049 39079 39081 39083 39121 39125 ...
## $ RUC_code : int 6 1 7 4 1 7 3 4 7 6 ...
## $ Pov_low_bo: num 13.9 11.5 7 12.9 15.2 14.3 14.7 8.5 12.7 8.1 ...
## $ Pov_Upper_: num 19.7 17.5 11 16.7 16.8 21.5 20.5 12.9 19.7 12.3 ...
## $ Pov_Percen: num 16.8 14.5 9 14.8 16 17.9 17.6 10.7 16.2 10.2 ...
## - attr(*, "data_types")= chr "N" "C" "C" "N" ...
str(slot(chi.poly,"bbox"))
## num [1:2, 1:2] -84.8 38.4 -80.5 42
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:2] "x" "y"
## ..$ : chr [1:2] "min" "max"
str(slot(chi.poly,"proj4string"))
## Formal class 'CRS' [package "sp"] with 1 slot
## ..@ projargs: chr NA
str(slot(chi.poly,"data"))
## 'data.frame': 88 obs. of 27 variables:
## $ OBJECTID_1: int 1 2 3 4 5 6 7 8 9 10 ...
## $ COUNTY_CD : Factor w/ 88 levels "ADA","ALL","ASD",..: 36 37 38 39 25 40 41 42 61 63 ...
## $ COUNTY_SEA: Factor w/ 88 levels "AKRON","ASHLAND",..: 33 42 51 58 20 35 70 53 10 61 ...
## $ ODOT_DISTR: int 9 10 11 3 6 9 11 5 10 1 ...
## $ FIPS_COUNT: Factor w/ 88 levels "39001","39003",..: 36 37 38 39 25 40 41 42 61 63 ...
## $ POP_2010 : int 43589 29380 42366 59626 1163414 33225 69709 60921 14645 19614 ...
## $ POP_2000 : int 40875 28241 38943 59487 1068978 32641 73894 54500 14058 20293 ...
## $ POP_1990 : int 35728 25533 32849 56240 961437 30230 80298 47473 11336 20488 ...
## $ STATE_PLAN: Factor w/ 2 levels "N","S": 2 2 1 1 2 2 1 1 2 1 ...
## $ ELEVATION_: int 1340 1220 1380 1200 1130 1040 1420 1420 1340 780 ...
## $ ELEVATION1: int 700 640 800 600 670 600 640 840 660 650 ...
## $ LAT_NORTH_: num 39.4 39.7 40.7 41.3 40.1 ...
## $ LAT_SOUTH_: num 39 39.4 40.4 41 39.8 ...
## $ LONG_EAST_: num -83.3 -82.2 -81.6 -82.3 -82.8 ...
## $ LONG_WEST_: num -83.9 -82.7 -82.2 -82.8 -83.3 ...
## $ AREA_SQMI : num 558 424 424 496 544 ...
## $ AREA_ID : int 53 54 55 56 143 57 58 59 75 76 ...
## $ SHAPE_STAr: num 2.41e+09 1.84e+09 1.91e+09 2.27e+09 2.40e+09 ...
## $ SHAPE_STLe: num 213158 215452 193045 197349 211754 ...
## $ NoOfUni : int 0 0 0 0 13 0 1 2 0 0 ...
## $ Shape_Leng: num 1.72 1.75 1.59 1.58 1.7 ...
## $ Shape_Area: num 0.151 0.115 0.117 0.138 0.149 ...
## $ FIPS : int 39071 39073 39075 39077 39049 39079 39081 39083 39121 39125 ...
## $ RUC_code : int 6 1 7 4 1 7 3 4 7 6 ...
## $ Pov_low_bo: num 13.9 11.5 7 12.9 15.2 14.3 14.7 8.5 12.7 8.1 ...
## $ Pov_Upper_: num 19.7 17.5 11 16.7 16.8 21.5 20.5 12.9 19.7 12.3 ...
## $ Pov_Percen: num 16.8 14.5 9 14.8 16 17.9 17.6 10.7 16.2 10.2 ...
## - attr(*, "data_types")= chr "N" "C" "C" "N" ...
summary(chi.poly@data$NoOfUni)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.500 1.295 1.000 14.000
plot(chi.poly)
leaflet(chi.poly) %>%
addPolygons(stroke = TRUE, fillOpacity = 0.5, smoothFactor = 0.5) %>%
addTiles() #adds a map tile, the default is OpenStreetMap
chi.poly@data$Pov_low_bo
## [1] 13.9 11.5 7.0 12.9 15.2 14.3 14.7 8.5 12.7 8.1 16.0 10.2 15.3 6.3 14.6
## [16] 7.5 9.6 5.3 8.1 14.9 5.0 13.3 11.2 9.1 13.6 16.9 18.3 25.7 11.5 7.3
## [31] 8.3 13.1 14.5 11.8 13.5 9.0 10.8 13.5 7.5 7.4 16.5 7.3 9.1 12.1 16.6
## [46] 7.5 16.8 14.0 5.0 16.4 11.8 17.0 8.8 12.4 7.3 7.6 12.9 6.7 15.2 5.4
## [61] 8.3 15.2 7.9 13.0 6.4 4.0 10.7 9.8 12.7 13.2 12.3 13.2 11.1 17.3 7.5
## [76] 8.8 9.4 7.6 4.1 10.3 7.2 12.9 10.5 15.7 3.8 9.7 9.6 9.4
require(RColorBrewer)
## Loading required package: RColorBrewer
qpal<-colorQuantile("OrRd", chi.poly@data$Pov_low_bo, n=4)
leaflet(chi.poly) %>%
addPolygons(stroke = TRUE, fillOpacity = .8, smoothFactor = 0.2, color = ~qpal(Pov_low_bo)
) %>%
addTiles()
chi.ols<-lm( NoOfUni~Pov_low_bo, data=chi.poly@data)
summary(chi.ols)
##
## Call:
## lm(formula = NoOfUni ~ Pov_low_bo, data = chi.poly@data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2727 -1.2258 -0.6412 0.1768 11.8641
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.20538 0.81222 -0.253 0.8010
## Pov_low_bo 0.13534 0.06905 1.960 0.0532 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.542 on 86 degrees of freedom
## Multiple R-squared: 0.04276, Adjusted R-squared: 0.03163
## F-statistic: 3.842 on 1 and 86 DF, p-value: 0.05323
list.queen<-poly2nb(chi.poly, queen=TRUE)
list.queen
## Neighbour list object:
## Number of regions: 88
## Number of nonzero links: 460
## Percentage nonzero weights: 5.940083
## Average number of links: 5.227273
W<-nb2listw(list.queen, style="W", zero.policy=TRUE)
W
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 88
## Number of nonzero links: 460
## Percentage nonzero weights: 5.940083
## Average number of links: 5.227273
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 88 7744 88 35.22186 356.7656
plot(W,coordinates(chi.poly))
coords<-coordinates(chi.poly)
W_dist<-dnearneigh(coords,0,1,longlat = FALSE)
moran.lm<-lm.morantest(chi.ols, W, alternative="two.sided")
print(moran.lm)
##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = NoOfUni ~ Pov_low_bo, data = chi.poly@data)
## weights: W
##
## Moran I statistic standard deviate = 2.0908, p-value = 0.03654
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I Expectation Variance
## 0.120686124 -0.014465249 0.004178391
chi.poly@data$chi.ols.res<-resid(chi.ols) #residuals ols
spplot(chi.poly,
"chi.ols.res",
at=seq(min(chi.poly@data$chi.ols.res,na.rm=TRUE),
max(chi.poly@data$chi.ols.res,na.rm=TRUE),
length=12),
col.regions=rev(brewer.pal(11,"RdBu")))
You can also embed plots, for example: